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We present a new scheme for the general computation of cosmic propagators that allow to interpo- 
late between standard perturbative results at low-fc and their expected large-fc resummed behavior. 
This scheme is applicable to any multi-point propagator and allows the matching of perturbative 
low-fc calculations to any number of loops to their large-fc behavior, and can potentially be applied 
in case of non-standard cosmological scenarios such as those with non-Gaussian initial conditions. 
The validity of our proposal is checked against previous prescriptions and measurements in numer- 
ical simulations showing a remarkably good agreement. Such a generic prescription for multi-point 
propagators provides the necessary building blocks for the computation of polyspectra in the context 
of the so-called F-expansion introduced by Bernardeau et al. (2008). As a concrete application we 
present a consistent calculation of the matter bispectrum at one-loop order. 



I. INTRODUCTION 



The large-scale structure of the universe that we observe today is thought to emerge from gravitational instabilities 
out of primordial metric perturbations, therefore precise observations of the large-scale structure of the local uni- 
verse can be used to put constraints on cosmological models. The connection between cosmological parameters and 
observations are however only straightforward when they correspond to the linear regime, i.e. when the observables 
can be computed as a linear transform of the primordial metric perturbations; this connection is less trivial when 
nonlinearitics, due in a large part to the gravitational dynamics itself, are present. Those nonlinearitics arise during 
the late stage of the gravitational instabilities. This is an epoch during which the universe can be safely assumed to 
be matter dominated (at least in the context of the standard model of cosmology where the dark energy component 
forms an homogeneous fluid.) 

According to the cosmological principle cosmic fields are statically homogeneous and isotropic. In the context 
we are interested in, two degrees of freedom are then relevant (see [l[ for details): the local density contrast, 5(p^) 
and the velocity divergence, 0(x) = 9j;Ui(x). It is then very fruitful to introduce the Fourier modes of the cosmic 
fields J(k) and 0(k), which evolve independently of one another in the linear regime. One can then introduce the 
doublet vE'a(k) = {i5(k), —6{k)/T-L], where "H is the conformal expansion rate with its time evolution described by the 
Friedmann equation, to write down the equations of motion in compact form and facilitate the implementation of 
resummation techniques. 

In this context, the goal of the theoretical and numerical calculations is a precise description, beyond the linear 
regime, of the statistical properties of VP^. We are particularly interested in the multi-component power spectra Pab{k) 
defined as, 

(*,(k)*,(k')) = <5D(k + k')Pab{k) (1) 

where the Latin indices a and b vary from 1 to 2, which are implicitly or explicitly measured in observations such as 
galaxy surveys, and higher-order spectra such as bispectra, 

(*,(ki)*b(k2)«'c(k3)) =(5D(kl+k2+k3)B„f„(ki,k2). (2) 

This problem is in general very complicated if one wants to solve it from first principles. It can be made slightly 
easier to address in the context of the standard cosmological model where the primordial metric perturbations are 
expected to follow Gaussian or near Gaussian statistics. In that case the primordial properties of the fields are entirely 
determined by the initial power spectra, Pa'b^'ik). The question is then to uncover the functional dependence oi Pab{k) 
as time evolves with P^"* (fc). 

In the last few years attempts have been made to present perturbative schemes in the context of the growth of 
structure. The standard perturbation theory is unambiguously defined but leads to uncontrollable results [l[. On the 
other hand alternative approaches have been proposed which produce more robust results, such as the Renormalizcd 
Perturbation Theory (RPT) in [2|, the Time Rcnormalization Group (TRG) approach in 0], the closure theory in 
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[Ij or with the help of perturbation theory expansion expressed in Lagrangian variables 0-0: etc. One of such 
approaches, advocated in [l^l, is based on what has been called the F-expansion. This method exploits the following 
relation, 

(3) 

where the standard Einstein convention (repeated indices are summed over) is used and where $a(fc) is the value of 
^a{k) at the initial time. This relation expresses the fact that ensemble average over primordial fluctuations can be 
reorganized in an alternative way to that in standard perturbation theory [l[ . It exhibits the multi-point propagators 
defined as the ensemble average of the infinitesimal variation of the cosmic fields with respect to the initial conditions. 
More precisely the multi-point propagators^ r^^''(qi, . . . , qp) arc defined as 

The relation ([3]) is valid for Gaussian initial conditions but can be extended for non-Gaussian initial conditions It 
clearly shows that the propagators are key ingredients for calculating nonlinear power spectra. They are the building 
blocks of the F— expansion approach and the focus of this paper. 

The second reason why these quantities appear to be important ingredients in perturbation theory calculations is 
that their asymptotic properties, e.g. how t hey behave for large wave-numbers, can be computed beyond perturbation 
expansions. This has been pioneered in 0, and recently revised in [l3j with the so-called eikonal approximation. 
More specifically it has been shown that in the high-fc limit the multi-point propagators are damped by a function that 
depends on the displacement field alone, irrespectively of the dynamics responsible of this displacement. Predictions 
on the behavior of those objects are then robust and can be computed in various approximations. 

This motivates their use as the building blocks of a perturbation theory scheme. To achieve this end, one must 
have a description of multi-point propagators at all scales, matching the perturbativc calculations at low-fc to the 
resummed asymptotic behavior at high-A:. As rcsummed propagators also contain information on loop contributions, 
both regimes can only be matched if a consistent interpolation schemes can be built. Furthermore, it is arguably this 
interpolation regime the most important for the prediction of loop corrections to (equal-time) correlation functions, 
as e.g. the power spectrum: while the low-fc limit can be safely computed using standard perturbation theory, the 
high-fc limit only adds a (time-dependent) phase-shift to Fourier modes and thus does not contribute to the power 
spectrum or bispectrum. 

This problem has been solved for the two-point propagator in [l^l with the help of an exponentiation scheme that 
interpolates between the one-loop results and its high-fc* behavior. However, this prescription is somewhat ad-hoc in 
that is specific to matching one-loop to resummed behavior of the two-point propagator but it is not clear a priori how 
it can be extended to incorporate higher-loop information at low-fc* or its generalization to multi-point propagators. 
The aim of this paper is to revisit this problem and propose a consistent solution which would be valid for any 
propagators and incorporate perturbativc information to any loop order. 

This paper is organized as follows. In the section II we recall the basic equations of motion, define our notation 
including the diagrammatic description and present the F-expansion of power spectra and bispectra in terms of 
multi-point propagators. In section III we review the expected properties of the propagators. In section IV the 
proposed interpolation scheme is presented in detail, while a comparison of our theoretical predictions for the tree- 
point propagator to numerical simulations is presented in section V. The implications of our results arc illustrated 
with a bispectrum computation in the section VI. Lastly section VII contains our conclusions. 

II. EQUATIONS OF MOTION AND THE F-EXPANSION 

A. The equations of motion 

We are interested here in the early stages of the development of cosmological instabilities in a cosmological dust 
fluid. In general the dynamical evolution of such a fluid can be described with the Vlasov equation for which one 



^ It is important to note that in this paper we call the F^*') functions the (p + l)-propagators - because it connects p -I- 1 lines - but 
alternative conventions can be found in the literature where F^^' is called the p-point propagator. 
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further assumes that multi-flow regions play a negligible role (see e.g. [13? 43 for recent discussion on going beyond 
this). In the single flow limit, the equations of motion then takes the form of a set of three coupled equations relating 
the local density contrast, the peculiar velocity field and the gravitational potential (see [![). 

At linear order these equations can easily be solved for an arbitrary background cosmology. One generically finds 
a growing solution and a decaying solution. Let us denote D+(t) the growing mode solution for the density contrast, 
with T the conformal time, and /+(t) its logarithmic derivative with respect to the scale factor so that, 

S{k,7^) = I?+(r?)<5o(k), 0{k,rj)/n = - f+{r^)D+{rj)So{k) (5) 

is the solution for the growing mode and similarly, 

5(k,r;) = D^{rj)So{k), 0(k,r,)/7^ = -/_(r,)i?„(7y)5o(k) (6) 

for the decaying mode. 

Following [l7|, the equations of motion describing a pressureless fluid in the one- flow limit can be written in a 
compact form with the use of the two component quantity 4'a(k, r), defined as 

^,(k,T)^[S{k,T), -j-^e{k,T)}, (7) 

where H = dlna/dr is the conformal expansion rate with a(T) the cosmological scale factor and where the index 
a = 1,2 selects the density or velocity components. Note that this definition of 'i'a is slightly different than the one 
used in the introduction and makes explicit use of the growing solution. The function /+ is unity for an Einstein-de 
Sitter background only. At this stage one can also remark that the choice of this basis is somehow arbitrary: we could 
have use any independent linear combinations of (5(k, r) and — 0(k, T)//+(r)7^ arc our choice of doublet fields. 

It is then convenient to rewrite the time dependence in terms of the growing solution and in the following wc will 
use the time variable rj defined as 

r, = log D+{t) (8) 

assuming the growth factor is set to unity at initial time. Then the fully nonlinear equations of motion in Fourier 
space (we henceforth use the convention that repeated Fourier arguments arc integrated over) read [l[ , 

d 

— ■fa{k,'q)+nab{v)'^b{k,7]) = 7abc(k,ki,k2) *b(ki,?7) *c(k2,?7), (9) 
orj 



where 

-1 

3 Sim 3 Sl,„ _ -1 
27f 



(10) 



and the symmetrized vertex matrix jabc describes the non linear interactions between different Fourier modes. Its 
components are given by 

n 1 1 ^ X n i i ^ I'^i +k2|^(ki ■ k2) 
7222(k,ki,k2) = dD(k-ki-k2) , 

7l2l(k,ki,k2) = fe(k-ki-k2) (11) 

7a6c(k, ko, kb) = 7ac6(k, kft, k^), and 7 = otherwise, where denotes the Dirac delta distribution. The matrix jabc 
is independent on time (and on the background evolution) and encodes all the non-linear couplings of the system. 
The formal integral solution to Eq. ^ is given by (see 0, [13, [l3 for a detailed derivation) 

*a(k,77) - 9ab{v) M^) + Cm 9ab{v,v') l^cdikM c{ki , v')^ d{k2, v') , (12) 

Jo 

where ^a{k) = ^a{k,ri = 0) denotes the initial conditions, set when the growth factor D+ = 1 and where gabiv) is 
the linear propagator, i.e. the Green's function of the linearized version of Eq. ([9]) and describes the standard linear 
evolution of the density and velocity fields away from their initial values. 

In the following calculations we will be using the value of the i^ab matrix to be that of the Einstein dc Sitter 
background thus assuming that = Qrn- Effectively it assumes that £>_ scales like D_^_^^^. This is known to be a 



4 



very good approximation even in the context of a ACDM universe. Within this approximation fiab becomes effectively 
time independent. It should be noticed that although the results presented below depend on this approximation, the 
whole construction is not based upon it. Calculations in an arbitrary background would simply make the whole 
presentation much more cumbersome, preventing the writing of explicit analytic forms. See Appendix A in Q for 
this. 

The ensemble average of any quantity can then be built out of the statistical properties of the initial fields. They 
are entirely defined from the initial power spectrum of density fluctuations Patik), 

($,(k)$b(k')) = <5D(k + k')Pab{k). (13) 

In what follows most of the calculations and applications will be made assuming initial conditions in the growing 
mode, for which $a(k) = (5o(k)ua with Ua = (1, 1), and therefore with PaT^'ik) = PQ{k)uaUb. 

The linear propagator gabiji) is one of the key ingredients and gives the variation of the mode amplitude as time 
evolves. The idea at the heart of the RPT approach is to generalize this operator beyond linear theory [l,[l3]- More 
specifically the quantity (5\l/a(k, 7y)/(5$b(k') expresses the way ^'a(k, ry) depends on ^^(k') as a function of time ry. 
This function however depends on the stochastic properties of the fields and one is led to define its ensemble average, 
GabiKVfiVt), as 

(^If^) - S^ik - k') Gabik, Vf.m). (14) 

where we have re- introduced the initial time rji. This quantity, known as the non-linear (two-point) propagator, 
depends on the initial fluctuations through the mode couplings. The ensemble average is made precisely over these 
modes. The Dirac-(5 function is due, as usual, to the homogeneity of the underlying statistical process. 

The expression for Gab can be computed order by order in perturbation theory. Such results can be obtained from 
a formal expansion of ^^(k, r;) with respect to the initial field, 

oo 

vE'a(k,r;) = VvE'(")(k,ry) (15) 



n=l 



with 



(k, r;) = (ki, . . . , k„; ry)$,, (ki) ...$,„ (k„) (16) 



where J^'"-* arc fully symmetric functions of the wave- vectors. Note that these functions have in general a non-trivial 
time dependence because they also include sub-leading terms in 77. Their fastest growing term is of course given by 
the well known G„} kernels in PT (assuming growing mode initial conditions), 

•^a6?b2...&„('^i'---'^"!'^)^''i •■•'"b" = 5D(k-ki...„) exp(™7) {i^„(ki,..,k„),G'„(ki,..,k„)} 

for a ~ 1, 2 (density or velocity divergence fields respectively). 

The concept of higher-order propagators is a natural extension of the non- linear propagator Gab- Such functions, 
that we denote F^^^^ (ki, . . . , kp), can be defined as, 

for second order (or three points), and for an arbitrary order they read, 

S <5i;:fl«7M > ^ - ''■■■■'>^S„.. (X.. . ^ . V, . (18, 

where ki p = ki -I- . . . -I- kp. They can be viewed as the building blocks of the theory. Note that for the purposes 
we consider here, we restricted our definition to derivatives with respect to the initial fields but a much more general 
description could be adopted. 

It is probably worth mentioning that F^^^^^ (ki, . . . , kp) are real functions which, for parity reasons, obey 

rS...., (-ki, . . . , -kp) = T^:l (ki, . . . , kp) . (19) 
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FIG. 1: Diagrammatic representation of the series expansion of ^a(k) up to fourth order in the initial conditions Time 
increases along each segment according to the arrow and each segment bears a factor gcdivf —Vi) if Vi is the initial time and rjf 
is the final time. At each initial point and each vertex point there is a sum over the component indices; a sum over the incoming 
wave modes is also implicit and, finally, the time coordinate of the vertex points is integrated from r]i = to the final time rjf 
according to the time ordering of each diagram. For instance, at fourth order there are two different possible topologies. 



C^,<'-'»p'(k,.i„n,) = 



!ag('l2-'l'2' 
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P,(q) 



1l 



gcb('i'rii) 



FIG. 2: Representation of the one-loop correction to the two-point propagator. The value of this diagram is obtained by the 
contraction of two incoming lines of 'I''''' multiplied by the initial power spectrum value. The expression of G'^^'°°'' is then 
obtained after integration over the internal indices c, . . . ,j, the momentum q and the times 771 and 772. 



B. Diagrammatic representations 

A detailed description of the procedure to draw diagrams and compute their values can be found in we can 
briefly summarize these rules here as follows. In Fig. [l]the open circles represent the initial conditions $b(k), where 
6 = 1 (5 = 2) corresponds to the density (velocity divergence) field, and the line emerging from it carries a wavenumber 
k. Lines are time-oriented (with time direction represented by an arrow) and have different indices at both ends, 
say a and b. Each line represents linear evolution described by the propagator gabivf ~ Vi) from time r]i to time rjf. 
Each nonlinear interaction between modes is represented by a vertex, which due to quadratic nonlinearities in the 
equations of motion is the convergence point of necessarily two incoming lines, with wavenumber say qi and q2, and 
one outgoing line with wavenumber q = qi + q2. Each vertex in a diagram then represents the matrix "fabc{<ii ^1,^2)- 
It is further understood in Fig. [1] that internal indices are summed over and interaction times are integrated over the 
full interval [0, 77/]. 

Loop diagrams appear once we calculate statistical averages such as correlators between fields. An example of 
such calculation (corresponding to the one-loop correction to the linear propagator) is presented in Fig. [2l where the 
average over initial conditions is encoded by the dependence on the initial power spectra P^"'' ((7), represented by the 
symbol (g). 

The previous construction can be extended to any number of loops. Note however that the explicit diagrams to 
be used depend on the statistical properties of the initial fields. For instance, for Gaussian initial conditions the 
two-loop diagrams contributing to Gab{k) are obtained from the contraction of 4 incoming lines in the expression of 
case the initial bispectrum is non- vanishing a non-zero contribution can be obtained from the contraction of 
3 incoming lines of 

These constructions can be pursued to higher-order propagators. A typical perturbation expansion of F'"^ is 
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r<"\k,pi,...,p^)= 




FIG. 3: Representation of the first two terms of tiie multi-point propagator F'"' in a perturbative expansion. F'"' represents 
the average value of the emerging nonlinear mode k given n initial modes in the linear regime. Here we show the first two 
contributions: tree-level and one-loop. Note that each object represents a collection of (topologically) different diagrams : each 
black dot represents a set of trees that connect respectively n + 1 lines for the first term, n + 2 for the second. 



presented in Figure [3l 



C. The F-expansion 

An important result rigorously shown in [Toj and recalled in the introduction, is that the series expansion of the 
power spectra can be rewritten in term of product of F'-P-' functions as, 

(*a(ki)^fc(k2)) =5D(ki+k2)^p! / d3qi...d3qp5D(ki -qi...p) 

x^it...a, (qi, ■ • ■ , qp) rit..fc, (-qi- • • ■ > -qp) Pa^tM) ■ ■ ■ pZ^'Mp)- (20) 

A further important property of this form is found when primordial metric perturbations are of only adiabatic origin. 
The initial power spectra then take the form P]Si^^ (k) = Ta{k)Tb{k)P'^^'^^-{k) where padiab.j-^-j jg primordial power 
spectrum of the adiabatic modes^. In this case, using the parity property from Eq. (|19p . we have 

(*4k)*fc(k')) = /d^qi..-d^qp[ri^al..a,(qi,---,qp)ra,(gi)T,^fe)" 

X [rit..^ (qi. ■ • ■,<ip)nAqimM] ^^^''^'-(91) • . ■p^'^'^'^-fe). (21) 

In this case, for a = the power spectrum is a sum of squares, that is, of manifestly positive terms. As for the 
RPT case, each of this term will correspond to a "bump" contributing to the final power spectra in a limited range 
of wavemodes. 

The resulting expressions depend on the primordial power spectrum and transfer functions. In this sense the result 
is a priori model dependent. Here, in the context of the F-expansion approach we concentrate on the late time 
behavior of the result thus keeping only the most growing terms. The linear transfer functions are then identical for 
the 2 components, Ta{k) = UaT{k) with Ua = (1, 1). We then simply have 

P^ {k) = Ua«6T2(fc)P-diab.(^) ^ UalibPoik) (22) 

and the expression pip can be rewritten in. 



(vl/,(k)vl/,(k')) /'d3qi...d3qpFiP)(qi,...,qp) F^^^ (qi, . . . , q^) Fo(gi) • • ■ Z'ol^p), (23) 

p P- J 

where Fi^-* (qi, . . . , qp) = Fia\...ap (qi, . . . , qp) Moi ■ • ■ Ua^- In the following we also assume that the structure of the 
vertices is that of an Einstein de Sitter universe. As mentioned before, it does not significantly restrict the validity 
range of these calculations. 



^ The latter that can be indifferently expressed in terms of the gauge-free potential perturbation or matter density perturbations. 
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FIG. 4: Representation of the resummation rule given by Eq. (|24|l . For Gaussian initial conditions, the bispectrum can be seen 
as a sum of product of r'*"' functions. 



This reconstruction scheme for the power spectrum in terms of F-functions can be extended to higher order spectra. 
For instance, generalizing Eq. (|20p the bispectrum can be formally rewritten by the resummation. 



(*„(ki)«'b(k2)*c(k2)) - J2 



r + s \ I s + i\ f t + r 
t 



r\s\t\ / d-^qi . . . d'^q^ d^q'^ . . . d^q^ d^q'/ . . . d^qj' 

\ /' / \ / \ t / 

r,s.t ^ ^ ^ ^ ^ 

x5D(ki - qi...r - qi...s) &(k2 + q'l...,, - q'/...*) &(k3 + q'/...^ + qi...r) 
xF('-+«) (qi, . . . , q., ql, . . . , q^) ri^+*' (-q^, . . . , -q',, q'/, . . . , qj') 

X F(*+'-) (-q'Z, . . . , -q;', -qi, . . . , -q,) Po(gi) ■ • • P^{qr) Po{q[) ■ ■ ■ Po{q's) Po{q'l) ■ • ■ Poiq't)- 

(24) 

This sum is diagrammatically represented in Fig. 2] We see that it runs over the number of lines that connect each 
side of the diagram (with the constraint that at most one of the indices r, s or t is zero, otherwise we would have 
a disconnected diagram). The leading order (tree) contribution is then obtained for r = s = 1, t = (plus cyclic 
permutations), up to one- loop corrections (in square brackets) we then have 

B{kiM,h) = 2F(2)(ki,k2)F(i)(fci)F(i)(fc2)Po(fci)i^o(fc2)+cyc. 

+ 8 / Ar(2)(ki-q,q)F(2)(k2+q,-q)F(2)(q-ki,-k2-q)Po(|ki-q|)Fo(|k2+q|)Po(g) 

(25) 



+ 6 j Ar(3)(_k3,-k2+q,-q)F(2)(k2-q,q)F(l)(k3)Po(|k2-q|)Po(<z)Po(fc3) + cyc. 
We will make use of this expression in section IVTl below. 



III. PROPERTIES OF THE MULTI-POINT PROPAGATORS 



A. The two-point propagator 

The large-fc behavior of the propagators can be addressed with the help of resummation techniques. This was 
pioneered in [l^l, taking advantage that for CDM spectra there is a characteristic scale set by the rms displacement 
(or velocity) field that sets the typical momenta inside loop diagrams, thus by large-fc we mean specifically k larger 
than this characteristic scale. This idea was put in a more general footing in [l3| where the concept of the cikonal 
approximation is introduced. In this context it is possible to compute the expression of the non-linear propagators 
taking into account the full resummed contributions of modes g <C A. 

The resulting expression of the propagator is then valid in the high k limit. More specifically one finds that 

n II. ^ . I \ ( '=^'^di>5pi.('^/''?*)\ /ORA 
Gab[k,'nf,r],) gabirif -i]^)cx.p\ 1 (26) 
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FIG. 5: In all the diagrams contributing to Gab{k), as the one depicted here, there is a line connecting directly the initial time 
to the final time. This is the principal line, drawn here with a straight thick solid line. The dominant loops contributing to the 
resummed propagator are those drawn by dashed lines, while the sub-dominant loops are those in dotted lines. 



where cTdispi. '7i) is the r.ms. of the one-point displacement field, d{T]f,r]i), from time rji to time rjf. More precisely 
the latter is given by 



assuming the velocity field is potential. The functional form (|26p is valid assuming the large scale displacement 
field obeys a Gaussian statistics. In that case the exponential damping is entirely determined by the variance of the 
displacement along one direction, 

'^Lpl.iVf,V^) = l{d^iv^■:V^))■ (28) 

In case the displacement is given by its linear expression and assuming it is dominated by the growing mode contri- 
bution one then has, 

<yLM^f:m) = (e"^ - e^'f J Poiq), (29) 

where Pa{q) is the initial linear power spectrum. This result was originally derived in [l^ from the explicit contribution 
of a large subset of diagrams - those that arc directly connected to the principal line (see Fig. [5]). The eikonal 
approximation shows that this result is very general. It is valid in particular irrespectively of the time dependence 
of the velocity field. As shown in [l3l |. this construction amounts to compute the displacement field from its linear 
expression. It is possible to include corrections to the displacement field statistical properties beyond linear theory. 
This was noticed in where 1-loop corrections to the variance of the displacement field are included in the calculation 
of the propagator damping function. It should be noted however that whenever the displacement field is not Gaussian 
distributed, the damping factor is not a function of its variance only. This can be naturally be taken into account in the 
eikonal approximation. For instance, the standard results can be extended to models with primordial non-Gaussian 
initial conditions for which one can recover both the resummation results and the T-expansion formulae (see [Tlj for 
details). In this case the exponential factor is replaced by 

exp ^ cxp(-g^(-^(-^^'g)-'-)^)-). (30) 

In this expression the variance of the ID displacement field along k has been replaced by the whole cumulant scries of 
the ID displacement field. Note this form can be extracted from another route, based on the use of Lagrangian space 
variables 0,3- this case, however, that it corresponds to the leading behavior in the higli-fc limit is ambiguous. 

On the other hand, as stressed in the introduction, the nonlinear expression of Gat can be approached with a 
perturbative series expansion. Formally one can write Gat as, 

Gab{k,rif,f]^) = gabivf ~ m) + ^G]^^"""^ {k^iqf ,7^1) + SGI^^"""^ {k,-qf + ... (31) 

where successive loop corrections are included. It is known that G^i^^°°'^ {k , i] f , rji) behaves like —k^a^^^p^ (rjf,rji)/2, 
etc., when k is large so that the perturbative expansion of (|26|) and agree when k is large. And this is expected 
to be true at all order in perturbation theory. This is actually the meaning one can give to the limit written in ()26p . 
Note that sub-leading terms of (^5)) are obviously expected to be present in the expression of Gab{k,r]f,r]i). If they 
appear within the exponential they would change the normalization factor. 

In all cases what we expect is that ([26| captures the large k behavior of the propagator whereas the expansion ((3T|) 
is expected to be a precise description of its low k behavior. One of the objectives of this paper is explore how these 
two limit behaviors can be matched together. 
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B. The RPT interpolation scheme for the two-point propagator 

The one- loop expressions of the two-point propagator have been expUcitly calculated in [l^ . We summarize them in 
this section, and their interpolation to the rcsummcd high-fc limit, as it will be useful to compare to our new proposal 
in section IV. 

The computation of the one-loop contribution involves in general different time-dependent functions. They arc all 
of the form a'^ where v is an integer or a half-integer. This is a consequence of the structure of the time dependence 
of the linear propagator and of the fact that we assume flab to be time independent. For each time dependence, 
each component 0^^^^°°^ has a specific k dependence that can be computed. But although there are 6 different a" 
functions that intervene in the exp ression of the one-loop diagram, the whole result can be collected into only 4 
different fc-dcpendcnt functions We recall here the explicit expressions of those results. One obtains, 



<5G}^'°°P(fc,a) 


= ^a{a) f{k) - 


^/3(a) i{k) - 
5 


h{k) - 


h ls{a) g{k), 
5 


5G}2^'°°P(fc,a) 


= ^a(a) f{k) - 


ll3{a) h{k) 4 
5 


- ^7(«) Kk) 


- ls{a) f{k), 


SGl^'°°P{k,a) 


3 

5 


h{a) t{k) + 
5 


3 

-7(a) i{k) - 


3 

-d{a) g{k), 
5 


dGl^'°°^ik,a) 


2 

= ga(a) gik) - 


lf3{a) h{k) - 


3 

-7(a) i{k) 4 


- l^ia) f{k), 



(32) 



with 

a{a) = a'^ 
and 

fik) = 
h{k) = 

m - 



7 2 



5 



—a — a - 
5 



^a-i/2, 7(«) = - a'/' + ^a-'^', '5(a) = 
5 5 5 5 



-a 

5 



-1/2 



-3/2 



1 



504fc3g5 
1 

168fcV 
1 

-1 
72fc3 g5 



6fc^g - 79/cV + 50g^fc^ - 21fcg^ + -(A;^ - q^ f{2k'^ + Iq^) In ^' 



^^0(9) d^q, 



^o(g) d^q, 



&k\ - 41/fc^<z3 ^ 2fc3g5 _ 3fcq7 + -(fc2 - q^f{2e + g^) In 

4 |fc + (7p 

6fc^g + k\^ + 9fcg^ + 7(fc2 - g2)2(2fc4 + r^k^^ + 3?*) In [^-=^1 Po(g) d^q, 

4 |K + 9rj 

Qk\ + 29k^q^ - 18fc3g5 + 27/cg^ -I- ^(fc^ - q^C^k^ + ^k^q^ + Sg"*) In i^— 4 

4 fc + a 



^o(g) d^q, 



2^0(9) 



We can notice that all these functions satisfy [l2[ 

/(fc), g{k), h{k), i{k) 

when k is large. The time dependent functions also obey remarkable properties, 

a(a) — /3(a) ~ a(a — 1)^, 



fc2 

2 



d3_q 
3^ 



(5(a) -7(a) = a-3/2(a_i)2^ 



(33) 



(34) 



(35) 
(36) 



so that, in the high k limit we indeed have 



5G^^'°°P(fc,r//,7?, 



:fc^crdispi.(?//>'7*) gab{Vf,Vi)- 



(37) 



One can also remark that a(a) is the most rapidly growing function and is therefore expected to dominate at late 
times. In this case only the functions f{k) and g{k) pla y a role in the expression of the propagators. Irrespectively 
of this limit, the proposed exponentiation scheme in |12| is the following. It is based on the exponentiation of terms 
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FIG. 6: Example of dominant loop contributions for the three-point propagator r^^^^(ki, k2, ks). The final expression is obtained 
by the sum of an infinite number of such loops and over all possible interaction times. 



either identified as the growing modes or the decaying modes, 





+ sg[\^ - 


-> Gii(fc, 


a) = 


3 

— ae 
5 


a,{a)f{k)-p,{a)t{k)) 


2 

+ 5^" 


-3/2g(5d(a)g(fc)-7i(a)h(fc))^ 


912 




^ Gi2(fc, 


a) ^ 


2 

- ae 
5 


a,{a)f{k)~l3,{a)h{k)) 


2 


-3/2g(5^(a)/(fc)-7d(a)h(fc))^ 


921 




^ G2l(fc, 


a) = 


3 

— ae 
5 


ag{a)g{k)+-ig{a)h(k)) 


3 

a 

5 


-3/2g(5d(a)g(fc)+/3d(a)i(fc))^ 


922 


+ SGi'^ - 


G22(fc, 


a) = 


2 

— ae 
5 


a,{a)g(k)-(3/2)^,(a)i 




3 a-3/2 g(5d(a)/(fc)-(2/3)/3d(a)h(fc)) 

5 



where we have redefined the a - ^ functions in Eq. p3p through aag{a) = Q;(a), aj3g{a) = a ^^^(3d{o.) = /3(a), 
a"fg{a) — a^'^/^7(i(a) = 7(a) and a~^^^ Sdia) = 5{a). At one-loop order these forms agree with the results (|32|) . 
They also agree with the limit form ([^5]) because of the properties p4ll36p . These forms also present a number of 
key properties: they are decaying functions of time and of k. This is ensured in particular by the fact that the 
terms under the exponential is always negative. Note that there are no free parameters in this construction: given 
an initial spectrum and cosmological parameters those form fully predict Gab{k,r]). These forms are the nonlinear 
propagator used in the RPT formalism. They have been successfully tested against N-body simulations so alternative 
interpolation schemes cannot significantly depart from them. 

It should be remarked however that those constructions do not necessarily represent the unique possible interpolation 
scheme. In particular if one allows the possibility of adding more than two exponentials, then one would obtain a 
whole set of alternative formulations. Before we move on to the formulation we propose in this paper let us first 
describe the multi-point propagator results. 



C. The multi-point propagators 



Let us continue with the diagrammatic approach, extended to multi-point propagators. The concept of principal 
line can be extended to the multi-point propagators. One can then define a "principal tree" which corresponds to 
the diagram when the propagator is taken at tree order (starting with order 4 there might be more than one possible 
tree). The diagrams contributing to the high-fc limit of the propagator are those that are directly connected to the 
principal tree [10|. 

For a given tree shape (q) (for instance, one of the diagrams of Fig. |6]), a careful resummation of all these diagrams 
gives the following result (for Gaussian initial conditions), 

T^(p-o) /, , / / \ / ^l...p'''displ.(^/' \ T^(p-g),trec/, , I I \ 

Tabi. ..6p(ki,---,kp, J?!, •••,?7p-i) ^ cxpl j Kbi.X (^^^■■■^K^Vf,Vt,Vi,---,Vp-i)- 

(39) 

where (q) labels a given topology and rjl are the time values at each vertex position. As this result is valid for any 
topology and any time, after proper summation we simply have, 

rip) n 1 \ ^ I ^^■■■p'^displ.i''^f'''^i)\ ^(p),troo/, , s 

Kb\...b^^^^■■■^K^Vf,V^) exp ^ j,^(ki,...,kp,7?/,77,). (40) 
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V : 



H8) 



FIG. 7: Large k, qi regularization of the one-loop diagram of the two-point propagator by higher order loops corrections in 
the high-fc limit. The main loop (in dashed line) is to be computed when q is in the UV domain. 



This result generalizes the one for the two-point propagator. Note that this result can also be derived in the context 
of the eikonal approximation showing that o'dispi.(j]fj'ni) does not need to be computed in the linear regime. 

Similarly to the two-point propagator, it is possible to obtain the low-fc behavior of the multi-point propagators by 
perturbative series expansion, 

^al...b, (kl, • • ■ , kp, 77/, 7?,;) = r^^^|'*;'°°(ki, . . . , kp, 7?/, 77,) + T^Pj;^'^"""^ (hi , . . . , kp, TJf, 77,) 

+rit'.t°'^(ki,---,kp,7?/,7;.) + ... (41) 

Note that in this case, even for the late-time behavior of the one-loop corrections, the relative sign between the tree 
term and the one-loop term is not fixed. It depends on both the geometry and the amplitude of the modes. 

Again, the aim of this paper is to propose an interpolation scheme between the known large and small scale 
asymptotic that fully respects these two limits. 



IV. PROPOSED INTERPOLATION SCHEME 



A. The case of the most growing mode 



The construction of our matching scheme is based on the analysis of the structure of the multi-loop terms corrections. 
To start with let us concentrate our presentation on the late time behavior of the propagator In this limit it is 
legitimate to assume that the initial fields are in growing mode only. We are then left with two independent quantities 
for the propagator, 

Ga+ = Gal + Ga2- (42) 

Up to one-loop order the late time expression of Ga+ is, 

Ga+ ^ e'^ + e^Vaik) (43) 

where fa{k) is either f{k), for a = 1, or g{k) for a = 2. The large k limit of fa{k) is fa{k) -k'^<7iispL/'2- The RPT 
expression of the propagator is then 

Ga+{k,7j)^e^cxp[e^ya{k)]. (44) 

The alternative prescription we propose is based on the following observation regarding the renormalization of the 
one-loop result: let us consider the diagram on Fig.[7]where the intermediate times r][ and 772 are fixed and the value 
of qi is also fixed such that its norm is large. The crucial observation we now make is that this object is then nothing 
but r'-'^^(qi, — qi, k, 77/, 77^), where the three incoming modes from the initial conditions correspond to the two dashed 
lines (joined at the initial power spectrum, i.e. the (x) symbol) and the rightmost /c-modc. But we now have a good 



In the context of the r-expansion approach it is in practice needless to consider sub-dominant terms and furthermore the extension of 
this construction to the general time does not introduce any new difficulty. 
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understanding of its renormalization properties, given by Eq. p9p . This equation tells us that we know how to resum 
all its loop corrections (some of which are represented by dotted loops in Fig. [7]) in the large qi and k limit. Because 
it corresponds to the effects of long wavemodes, let us call this the infra-red (IR) correction. The expression we find 
is an application of the general result for multi-point propagator and it leads to a simple exp(— fc^cr^j^pj /2) factor. 
It is important to note that it is independent of qi (and intermediate times rji and 7/2)- We still have to perform 
the integral 771 and 772 and then over qi . For the latter however we have to bear in mind that it cannot run over all 
possible values: it has to avoid its IR part. We can then split this integral into 2 parts, one IR, for which this result 
is not valid, and one UV for which it applies. Then the value of that set of diagrams would be 

<5G.+ = e^'^f^'^Hk) + e'Vr\k) exp (-fc^aLp,./2) ■ (45) 

We then note that the first term is simply the first (non trivial) term of the usual IR resummation of diagrams, 
exp(— A:^crj;gp[ /2). We are then let to simply set 

/i'^Hfc) = -^fc^-Lpi.. (46) 

/r^^^fc) = fa{k) ~ fi'^\k). (47) 

This identification leads to the following form, 

G-f = eJ' + SGa+ = e" (1 + e'\Uk) + fcVLp,./2) exp (-fcVLpi./2) (48) 

for a "regularized" propagator which compared to the expression (|44[) amounts to replacing 

exp (e^Valfc) + fcVLpi./2) ^ (1 + e^'^faik) + fcVLpi./2). (49) 

Note that these quantities are both finite at large k. In Fig. O we compare these two prescriptions for the density 
and velocity fields at z = for a A-CDM cosmology. They are virtually indistinguishable when one considers the 
propagator shape. They significantly depart from one another only for fc ~ 1 Mpc, as shown on the right panel, 
that is at scales where the exponential damping is already extremely strong. 




Jh , , , , , , , , I ^ L , . I . . . . , . I .... I ^^ 

0.01 0.02 0.05 0.10 0.20 0.50 1.00 0.01 0.02 0.05 0.10 0.20 0.50 1.00 2.00 

k k 



FIG. 8: Comparison of prescriptions for exponentiation schemes for the density (black upper lines) and velocity (brown lower 
lines) fields. The thin solid lines correspond to the RPT prescription, the thick dashed lines to this work. The left panel shows 
the propagators, the right corresponds to the propagator when they are divided by the global damping factor exp( — fc^adispi./2). 
The computations are made for z = and for a A-CDM cosmology (see Section |V] for details). 



B. The general prescription 

Since the resummation properties of the F functions preserve the topology and the intermediate time structure, the 
whole procedure applies to the full time dependence of the one-loop term. The integral splitting can then be done 
more generally and one gets. 



displ 



exp 



Mispl 



(50) 
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The relation with the RPT proposed forms given by Eqs. (|38|) is here more comphcated. It is clear however that both 
propositions agree at the one-loop correction level and both propositions exhibit the same high-A: behavior. 

One important aspect of this construction, which will be exploited in the following, is that it can obviously be 
extended to multipoint propagators, 



+ i^k^(^dispiiVf,m)K'bi...b^0^i,- • ■ , kp, 77/, 7?,) 



exp 



(51) 



where fc = |ki + • ■ • + kp|. By construction this form is such that it has both the correct one-loop correction and 
the correct large-fc behavior. One can also remark that, unlike the RPT interpolation proposal where the matching 
is done in the "basis" given by density and velocity fields, the resulting form here is independent on the basis chosen 
to represent the cosmic fluids. This may not be an important difference for the CDM-only case but could be of 
importance when one wishes to extend the field content to other degrees of freedom (e.g. the case where one has in 
addition baryons, massive neutrinos, extra fields in the context of modified gravity theories, etc.) 

Finally, another important advantage of this construction is that it can be pursued to any order in loop corrections 
in a rather straightforward way. 



" abi ...bp 



-ptrcc 
^ abi.. 



6T 



1 — loop 
abi ■■■^p 



1 



7,2 2 ptrcc 
^ "^displ.^ abi.. 



■ srx:z + c.t. 



exp 



Mispl. 



(52) 



where c.t. is a counter-term such that the 2-loop expression of the reg. expression is exact. 



c.t. 



displ. 



ntlCC 

■ abi ■ ■ 



displ. 



rpl-loop 
abi. ..bp 



(53) 



Let us now illustrate this with some examples. The resulting interpolation scheme for T^"^^ for specific geometries 
is shown in Fig. [9] for its late time behavior, showing that the interpolation is rather smooth. In particular it can 
handle the fact that tree-order and one-loop corrections have different signs. This is the case illustrated in the lower 
right panel in Fig. |9l 



C. The case of Non-Gaussian initial conditions 



The case of PNGs can similarly be taken into account. In this case the damping factor is changed in order to take 
into account the higher order cumulants of the ID displacement field as given in Eq. (|30|) . Still it is possible to apply 
the same regularization scheme except that the counter terms have to be recomputed. 

Novel two-loop order terms, depicted on Fig. 1101 are due to the primordial bispectrum. In the eikonal approximation 
(e.g. when the vertex values are computed for qi <?C fc), these diagrams vanish however. They therefore do not have 
counter terms in the regularization scheme we propose. Differences in the counter terms appear only at the three- 
loop order. It corresponds to the fact that the damping factor given by Eq. pO[) is not sensitive to the primordial 
bispectrum, but it is to the primordial trispectrum. 



V. COMPARISON WITH NUMERICAL SIMULATIONS 



As the prescription we are advocating here does not give significant differences for the two-point propagator (already 
studied in [13) we focus our analysis in the three-point propagator. 

As shown in the previous paragraph the prescriptions in Eqs. (|5T| [52]) give non-trivial behaviors for the three-point 
propagator. These prescriptions can be compared against measurements in numerical simulations provided one can 
measure cross-bispectra between initial conditions (5o(k) and the final density fields J(k, 77). The estimator for the 
three-point propagator was introduced in [loj ; 



r 



^"\h,k„k,) ^ \ 1 E E^(k^^) X '5o(-k05o(~k,), (54) 



' ™ 2Po{k,)Po{k,) N ^^^^ ^ 
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FIG. 9: The shape of the r'^-* propagator for "equilateral" (fci = ^2 = fca = fc, top left), "colinear" (fci = k2 — k/2 ; kz = k, top 
right), "anti-colinear" {k\ = 2k ; k2 = ka = k, bottom left) and "squeezed" (fci = k2 = ik ; kz — k, bottom right). The long 
dashed (blue) lines correspond to the simple exponential cut-off obtained in the high-k limit, Eq. (|40p . the short dashed (violet) 
lines are the tree plus one-loop order results and the solid (brown) line correspond to the interpolation scheme proposed this 
paper. These comparisons are made for a power spectrum corresponding to A— CDM model at 2: = 0. 



^ bispecti'um/1 , k ^ , k ^ , k , 

Gab w= — I 1 — — + — ^ < / + — ^-7 — 

-^-(^^^ (XX__>' 

Bp(q,,q2,q3) 

FIG. 10: Dominant complementary diagrams contributing to the one-point propagator in case of PNGs. They make intervene 
the primordial bispectrum Babc{<li, 92, gs). For models of cosmological interest they are of the order of the 2-loop contributions. 
When those diagrams are computed in the eikonal approximation, that is when the vertices are approximated by their high k 
limit they all vanish for parity reasons. They do not so in general. 



where the sum runs over Fourier modes in the |ki| bin, kj in the |k2| bin and k; in a bin jkaj such that |ki — k2| < 
1^3 1 < ki -|- k2, and N is the number of terms in the triple sum. Equation (|54p is valid for initial conditions set in 
growing mode for which the only measurable quantity is the contraction T abib2'u-biUb2 with u = (1, 1) (and a = 1 = 6 
in our case). In addition it assumes Gaussian initial conditions, see [ll| otherwise. We note that a similar expression 
holds for the velocity divergence propagator (i.e. a = 2), but its study is beyond the scope of this paper. 

We used Eq. ([51)1 to measure the three-point propagator in a set of 50 N-body simulations, each containing Npar = 
640'^ particles within a comoving box-size of side Ltox = 1280/i^^Mpc. The total comoving volume of our simulations 
is approximately 105 {h~^Gpc)^ . Cosmological parameters were chosen as flm = 0.27, JIa = 0.73, fit = 0.046 and 
h = 0.72, together with scalar spectral index = 1 and normalization dg = 0.9. The simulations were run using 
Gadget2 [2^ with initial conditions set at Zi ~ 49 using 2nd order Lagrangian Perturbation Theory (2LPT) [l3, Hlj . 

Before comparing theoretical predictions and numerical results it is important to account for binning effects since 
correlations of modes within bins implicitly change the shape of the predicted high-order propagators. Hence we will 
proceed by computing the predictions for binned modes. This is easily accomplished as follows. Let us denote by an 
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0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 



k [h-^ Mpc] k [h-^ Mpc] 




0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 

k [h-^ Mpc] k [h-^ Mpc] 

FIG. 11: Comparison of the proposed interpolation scheme against measurements of the two-point density propagator in 
numerical simulations at z = 0. Solid line is the prediction from Eq. (|51|l including the only most growing mode into p(2)troo 
and 5p(2)i-ioop g^fj-gj. jg corrected for binning as discussed in the text; the dashed line is when the binning is not taken into 
account and the dotted line is the prediction when the one- loop correction is not taken into account. Symbols with error 
bars are the corresponding measurements in 50 runs of Lbox = 1280/i~^ Mpc each. Different panels corresponds to different 
wave- vector configurations (or "triangle shapes") as detailed in the y-axis labels (and main text). The results are plotted as a 
function of k — ks. 



overline the bin average, e.g. 

r^'^(fcl,fc2,fc3) = -— ^-^=— / d^qi / d3q2 / d3q3P(gi)P(<Z2)r(2)(qi,g2,<Z3)'5D(q3-qi-q2) (55) 
NhinP[kl)P{k2) JBi Jb2 JBs 

where q^ is in bin Bi = {ki — 5k/2, ki + Sk/2} and A'bin is the normalization 

^bin- / d^qi / d^qa / d3q3 (5D(q3 - qi - q2), (56) 

JBi JB2 JBa 

and 

P{h) = / d3q,P(g,)/ / d^q,. (57) 
Jb, JBi 

Then writing the Dirac ^-function as (5D(k) = J ^^^"3 exp(— ik.u) we finally have, 

T^^\ki,k2,k^) ^ I u^du [ ql io{qiu)dqi [ qi jo{q2u)dq2 [ 93 jo(93")'^93 r(^'(gi, 52, 93) (58) 

^^bin Jo JBi JB2 JB3 
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with iVbin « 8Tr'^kik2h{Skf. This expression, where F^^^ is the model prediction from Eq. (|5ip . is the one we use 
to compare with measurements obtained in numerical simulations. Notice that we are now using three wave-modes 
as arguments with ka = ki + k2 being the "outgoing" momenta. More precisely, in our computations we choose the 
central bin value to be fc^ = 4i x 27r/Lbox and the bin width to be fixed and given by (5fc = 4 x 27r/Lbox- Specific 
geometries, i.e. ratio of wave modes, are then defined from the central values of the bins. 

Those comparisons are shown in Fig. [TT] for four different shapes, "almost squeezed" (fci = fc/4 ; fc2 = ^3 = fc), 
"elongated" (fci = fc/2 ; k2 — fcs = k), "equilateral" (where ki = k2 = k^ = k) and "colinear" (fci = ^2 = fc/2 ; k^ = k). 
(Note that these configurations are not the same as in Fig. O This is due to the fact that the "anti-colinear" and 
"almost squeezed" configurations presented there exhibited too large signal-to- noise ratios). Remarkably all measured 
configurations show a very good agreement between the numerical results (points with error bars) and the proposed 
prescription (solid lines). The dashed lines show the prediction before the binning corrections. The latter can be 
quite large specifically at large scale (as expected) and furthermore in some configurations one-loop term corrections 
significantly improves the predictions when compared to the numerical results. One can also note that the size of 
the error bars change from configuration to configuration. This is due to the fact that, for a given fcs the number of 
available modes in all three bins vary strongly. 

Overall, these results indicate that the proposed interpolation scheme works remarkably well when compared to 
measurements in simulations in an extended fc-range, which in turns is a very important step towards the accurate 
modeling of polyspectra. 

VI. APPLICATION: CALCULATING THE BISPECTRUM 

We are interested in the computation of _Ba6c(ki, k2, ka), the bispectrum of the components a, &, c of the cosmic 
fluids, 

(*a(kl,?7)*;,(k2,77)^'c(k3,77)) =fc(ki+k2+k3)Babc(kl,k2,k3,77) (59) 

Such bispectra can be computed from a resummation of product of F functions. This is an extension of and this 
is illustrated in Fig. [T^l In particular if one wants to incorporate all one-loop effects one should include three types of 
diagrams: the first one involves the one-loop correction of the propagators, the second, third and fourth correspond 
to intrinsic one-loop contribution to the bispectra. 




FIG. 12: The contributions to the bispectrum up to one loop order. In the first diagram the expressions of r(2' and r<i' should 
incorporate their own one loop correction ; they are computed using formulae (|5ip at one loop order. In the other 3 diagrams, 
the propagators correspond to their regularized tree order expressions; they are computed using formulae (|40p . There is thus 
no need at this order to include their one-loop correction, in particular to F'^-*. 



More precisely, for growing mode initial conditions, the bispectrum takes the form, 

BaU^iMM.v) = 2Fi'|^-^°°P(k2,k3,r;)G,^;'°°P(k2,7y)G^;'°°P(k3,,y)Po(fc2)Po(fc3) +sym. (2 terms) 
+ 8fc(ki + q2 - q3)(5D(k2 + qs - qi) ^(ks + qi + q2) ^^0(171) ^0(92) ^0(93) 

T^(2)trcc/ \Tn(2)trcc/ \-|^(2)trcc/ n , /o j_ \ 

X F;^+ (q3, -q2, 77) F^^'_^ (qi, -qg, 77) T\^^ (qa, -qi, ??) + sym. (2 terms) 

-f 12 Fi^|^'^<=(-k2 + q, -q, -k3, 77) Tfl';''\\^2 - q, q, v) (^3, v) 

X Po(|k2 - q|) Pa{q) Poiks) + sym. (2 terms) (60) 

where all the propagators are computed using their regularized form, whether it is at tree order as in (|40p or including 
one-loop correction as in (|5ip . Integrals over the wavevectors qi, q2 and q3 when present are implicit, and the 
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symmetric terms are obtained by a simultaneous permutation of the indices abc and the wave modes ki , k2 and ka . 
Note that in this formulation the T^p^ functions are assumed to be symmetric functions of their arguments. Because 
of that the last 2 diagrams that appear in Fig. [12] are automatically included. 

(2) 

For its practical evaluation the non trivial part is the one-loop expression of the r^^_|_ propagators. Its properties 
were discussed in [l^. We give in the appendix their actual values for both the density and the velocity fields. They 
are the crucial ingredient to use to compute the bispectrum at one-loop order. 
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FIG. 13: The shape of the density bispectrum, _Biii(fci, ^2, fcs), as a function of scale for different configurations of the 
wavevectors. The (gray) dotted line is the standard tree order prediction for the bispectrum and the (gray) dot-dashed the 
standard one-loop result; the (green) short dashed line corresponds to the first contribution in the diagrammatic expansion of 
Fig. 1121 the (red) long dashed line to the 2 others and the solid line the resulting reconstructed bispectrum. The first panel 
shows B{k, k, k), e.g. equilateral configuration. The second corresponds to a lightly squeezed configuration, k'^ B{k/2, k, k), 
and the last to a flatten configuration, k'^B{2k, k, k). The predictions are made for z = 0.5. 




e 

FIG. 14: The reduced bispectrum as a function of the relative angle between the wave- vectors 9 for fci = 0.1/i/Mpc and 
k2 = 0.2/i/Mpc. Lines follow the same convention as in Fig. 1131 



We present in Figs. [T51 and [Til the resulting shape of the density bispectrum for a standard ACDM model &i z ^ 0.5. 
Figure [T5I shows the scale dependence of the bispectrum (multiplied by to make it less scale dependent). It makes 
clear that the contribution of the first diagram and that of the 3 others correspond to different scales, each producing 
one bump at different scales. This is reminiscent to what the RPT calculations give for the power spectrum 

In turn. Fig 1141 explores the resulting angular dependence of the reduced bispectrum. Here we plot 

Q(fci,fc2,fc3)= „ n \v II, \ , V (1 \u (1 \ , u II \v II x ^iii(ki,k2,k3) (61) 

PQ{ki)Po(k2) + Pa{k2)P^{ki) + Po(fc3)^^(Kl) 

where the expression of the power spectrum is kept at linear order. The result is expressed for fixed values of 
ki = 0.1 h/Mpc and ^2 = 0.2 h/Mpc and as a function of their relative angle, 6, so that fc| = fcj -I- fc| -I- 2fcifc2 cos6'. 
The plot compares the naked "tree" order result (dotted line) with our prediction or the perturbative calculations 
(same convention as in Fig. I13p . Detailed comparisons of such predictions with A'^-body results is left for future 
studies. 
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VII. CONCLUSIONS 

111 this paper we propose a systematic interpolation sclieme that aims at describing the (multi-point) propagators 
in such a way that their expressions interpolate between their perturbation theory forms - to any loop order - and 
their large-fc behavior obtained from non-perturbative re-summations. This scheme is based on a separation of scales 
between the long wave-modes, whose effects are fully resumed and lead to the large-A: behavior, and the short wave- 
modes that are fully taken into account in the perturbation theory treatment. Our prescription is different from the 
exponentiation scheme proposed in [l^ but departs from it only very weakly. Our new prescription however is very 
general and can be used for any loop order calculations and for any propagator. Furthermore, this construction is 
totally unambiguous. We note that it can even be used in the context of primordial non-Gaussianities although new 
terms arise only at two loop order. 

As the construction proposed here has a very general range of application, it should in principle be tested for a large 
variety of quantities, from two-point and multi-point propagators of the density field to the ones for the velocity field. 
We proposed here some comparisons with N-body results for the quantities that are of most interest for the use of the 
multi-point propagators in the context of the F-expansion applied to the calculation of the density power spectrum. 
Because of the property we mentioned in the previous paragraph, this prescription for the two-point propagator is 
found to give very accurate results for the two-point propagator. We leave for further studies the impact of two-loop 
effects. In this work we further check the validity of our prescription against numerical results for the one-loop level of 
the density three-point propagator. These comparisons are presented in Fig. 1111 We found that it gives a satisfactory 
form for a large range of configurations, i.e. interpolating the low-A; one-loop result with the high-fc exponential decay, 
even when the signs of tree order and one-loop forms differ. 

In the context of the F-cxpansion approach this construction therefore provides us with the necessary recipes for 
constructing poly-spectra incorporating any order of perturbation theory results. In coming papers we will explicitly 
compare predictions of the F-expansion for the power spectra with numerical simulations when higher order FT loop 
corrections are included. These prescriptions provide a good opportunity for giving the explicit form of the bispectra in 
the context of the F-expansion. The explicit mathematical forms are given in Section IVll when bispectra are computed 
up to one-loop order. Such expressions make use in particular of the three-point propagators at one-loop order whose 
explicit forms are given in the appendix. We found in particular that the bispectra terms can be separated in a tree 
order contribution and coupling terms contributions that contribute at different scales, i.e. in subsequent bumps, in 
a similar way to the power spectrum. Comparison of the proposed forms for the bispectrum with iV-body results is 
however left for further studies. 

We finally note that such a construction is restricted here to the case of a single pressure-less fluid. Whether it 
could be used in cases the fluid content of the universe is rich er ( with extra degrees of freedom carried by baryons, see 
[13I. [2pl . massive neutrinos, dark energy components, see [H, [13; etc.) is still largely an open issue. Results obtained 
in [1^ however suggest that the effects of adiabatic long wavelength modes can be resummed and therefore could be 
incorporated in a large variety of cases. 
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Appendix A: Explicit expression of the one-loop F'^^ propagators 

We give here the explicit expressions of the late time component (most growing term) of the 1-loop expression of the 
three-point propagator (for growing mode initial conditions). The final result is obtained by the explicit integration 
over the wave mode q depending on the linear power spectrum Po{q), 

r(^^^-'°°P(fci,fc2,fc3) = W J q^dqT^^l'-'°°''ih,k,,k,,q)Po{q) (Al) 

where the result is expressed in terms of the 3 norms, fci, fc2, k^ so that fc| ~ fcj + fc| + 2ki.k2. Note that these three 
norms obey the triangular inequality. 

Note that r|^^^ '°°''(fci, k2, k^, q) obeys the following asymptotic behaviors, 

r(^^r'°°'(fci'fc2,fc3,g) ^ -^r(^'r(A;i,fc2,fc3), (A2) 

and q ^ 0. We also have 

rf^^i-'°°p(ki,k2) ^ (A3) 

r(^ir'°°P(ki,k2) ^ ^9ik2), (A4) 

when ki — >■ and where the functions f{k) and g{k) are defined in Eqs. ([55)1 . These 2 asymptotic regimes can be 
used as internal checks. 

The function T^^^ ^°°^(ki, k2,k^, q) has been obtained after integration of the angular variable in the loop expres- 
sions. The result can be expressed in terms of the functions 



(fcl - qf 



^(^i) = log (A5) 



^ . -^^3.- (fej - - ^hqVi-k-t - k-i + kj) q^ + k-tk-i + q^ 

\ -Akiq^ - 2 (fc2 - g2) (^2 _ ^2) + 4yfc3gV(-fc^ - k'i + ki) g2 + kfk'i + q^ J 
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and it reads, 



t(2)1-1oop 
1++ 



{ki,k2,k3,q) 



k2 {kl 



2fc| 



29568 fc^gV?^ (k'i ~ki+ g2) + kf {k'i - 



{k 



17248 klklq^^q^ {-kl + kl+ g2) + kl {k'i 



W{ki,k3,k2) +sym. (ki ^ 



W(fci,fc2,fc3) 



165 {kl - kj) ~ kl° {582kl 



1655808 kfklq^ 
+kl (-lllSfc^g'' + 179k^q^ - 3kl (211fc. 
+kf (I787fc2g2 + kl (2120A:| + 6891?^) - 
+k^q^ (405fc^q'' + kj (7269g2 - I677fc|) 
-3 {kl - kl) klq^ {kl {285q^ - 473/.|) - 

(7fc|g2 + 2fc| - 9q4 



79U-3 + 1947g2 



2155g^ 



kl (758fc2g2 



390A:? - 7284g*) + 492/c° + 151fcS + 219g^ 



51/c| 



- 321kjq^ 
17 (gfc^g^ 



601fc;^ 4 



1689g*) 
-4376A:|g2 
) + 558fc|) - 



- 391fc| 
39/cf] 



+ 471g^) + 1017A:^ + 269fc^) 
C{ki) +syni. (ki ^ k2) 



236544 klklklq'^ 

15 (fcj + fc: 



-12(fc?-fc22)%«-3(fc?-A;|)'(fc2 



kl) q' 



-2kf 



2) ''^f 



16(fc2 



^2^ 



-fc|H7(fc2 + fc2)g4 



lkt)q^ 



{kl 



2q 

kl) {kf 



(46 {kl 
8 (7A:^ - 17 klkl 

+kl {kl + kl) q^ + {52kt - 96klkl + 52fc|) q'^ + 2 {kl - kl)" {kl + 
1 



lOfc^fcJ 
kl) q 



C{k^) 



{mklklkl^ 

iOQklkl {kl 



kl)q' 



6209280 q4fc8fcffc| 
+ [-2475 {kf 

-68910fc?fc^(?2 + 5415fc?fcf {kl + kl)] k 
+ [135 (55fc^ - blklkf - h\k\kl + 55/c^) 



10 
3 



9090/cffc^ (fc2 + fc2) q4 
45fc?fc| (fc? + (663fc]^ - Sllfc^fcJ + 663/c^) q^ 



+ (38925fc|fc? + 9532/c^fc? + 38925fc^fct) q"- - 120370fc?fc^ {kl + fc|) 
+15fc?fcf (377fcf + 940fc|fcJ + 377/c|)] kl 

[-135 (fc? + fc^) (55fci* - 201fc|fc^ + 156fc|fcf - 201fc^fc? + 55fc^) 



+30fcjfc^ (1959/c^ - 3283fc|fc^ - \9Qlk\k\ - 3283fc^fc? + 1959fcf) (t 
-2k\k\ {kl + (53220fc]^ - 178337fc|fcJ + 53220A:|) q"- 

+10/c?fcf (22261fc^ - 45566fc|fcJ + 22261fc|) g^ - 15A:^fc^ {kl + /c|) (805fct - 2306A:|fc? 
[45 (55^1^2 - 28hklk^^ + 157A:^fc^ + 230/fc^fc? + 157A:^fc;' - 285A:^°fc^ + ^^kf) q^ 
-Ihklkl {kl + kl) (l949fc^ - 10643A:|A:^ + 20573fc|fc]' - 10643A:^fc^ + 1949fc^) q^ 
+3fe]'fc| (19475A:^ - 139742A:|A:? + 246974fc^fc;' - 139742/c^A:J + 19475/c^) q^ 
-iQklkl {kl + /cf) (911fct - IhlQklkl + 911/0^) q^ + 45A:?A:^ (fc? - kl) ^ (l3A:f - 4fc^A:J 
7560fc^fc^ {kl + /.|) ^8 + \2mklkl {A?,k1 - 80fc2fc2 + Aikt) q^ + 3150fc?fc^ {kl ' ^ 

11340g*^fcffc^ {kl - fc^)^ - 2835g^fc^fc^ {kl - fc^)^ (fc? 



805fc|)] fc. 



13^2^)] fc; 



(A7) 



It can be noted that the final expression is symmetric in fci and fc2. 
constructed, 



The velocity component can be similarly 





^ [55 (fc| - kj) _ kl° (l94/s| + 1921fc^ + 6i9q^) + 13kl 



+kf {25klq^ - n59k^q'^ + k^ {SOlkj - 2155^^) + k^ (-3414fc|g2 + 306/c| - 2428g**) + IdAk^ + 329fcf + 73q^ 
+kf (2913fc^g2 + kj (ll76A:| + 2297g2) + 17kl + 1579k^ + 563q'^) 

+k\q^ (l35A:3g^ + k^ (321/c| + 2423?^) + 333A:;^g^ + kl {-IMklq^ - 1307^3 + 1579"*) + 339/c^ + 647fc^) 
-3 {kl - kl) klq^ (kl (I2lkl + 95g2) - ^Iklq^ + 186A:| - 307fc3)] C (fci) + sym. (ki o ks) 



+ (-825 (kl + kl) q^ ~ UUOklkl (kf + k^) q^ + 32290ktk^ (kj + kl) q^ 

-133850fc^/c^g2 ^ u^sbkfkl (kf + kl)) kf 
+ (45 (55fc^ - blklkl - 51fc|fc2 + bbkl) q^ + 15/fc?fc| (/cj + kl) (lOQ9kl - A2lklkl + 1009/s|) 

-fcj'fc^ (50385A:^ + 92816A:^fc2 + 50385fc^) - 81190fc^fc^ (kl + kl) q^ 

+l?,klkl (683fcj' + 1252klkl + 683k^)) k^ 
+ (-45 (kl + kl) (55kl - 201klkl + 156k*kf - 201fc^fc2 + 55fc^) q^ 

+30klkl (235kl - 757klkf - 137k^kf - 7h7klk\ + 235fc^) q^ 

-^k\k\ (k\ + kl^ (3A5,k\ + 628A:^A:^ + 345A:|) q^ 

+Wik\kl (24507fcj' - 50506fc^/fc2 + 24507A:^) - \hk\kl (k\ + fc^) (\977k\ - mWiklk\ + 1977fc^)) kl 
+ (15 (55fcf - 285A:^fcf + 157fc^A:^ + 230fc^fc^ + 157fc§fc;' - 2%hk}^k\ + hhk)^) q^ 
-5klkl (kl + kl) (I949fc^ - 10643fc^fc^ + 193Alk^kf - 10643A:^fc^ + 1949fc^) g*^ 
+kfk^ (19475A:^ - 93542A;^fc^ + 160734fc^A:^ - 93542/c^fc^ + 19475fc^) q^ 

-90klkl (kf + kl) (67kl - lOGklkl + 67k^) q^ + Wkfk^ (kl - kl) ^ (l3k^ - Aklkf + 13/s|)) kj 
+ (2520klkl (kl + kl) + UGOk^k^ (7kf - 12klkl + 7k^) q^ - 1260/c^fc^ (kl - kl) ^ (kl + kl) g^) kl 
-3780g^fc^/c^ (kl - kl)^ - 945g^fcjfc^ (kl - kl) ^ (kl + kl) } . 



+ 





"-12 (kl - klfq^ - 3 (kl - klf (kl + kl) g4 



+ 



